A Method for Characterizing Hydrogeologic Heterogeneity Using Lithologic Data Gregory P. Each1, L. Larry

Large-scale (> 1 m) variability in hydraulic conductivity is usually the main influence on field-scale groundwater flow patterns and dispersive transport. Incorporating realis tic hydraulic conductivity heterogeneity into flow and transport models is paramount to accurate simulations, particularly for contaminant migration. Sediment lithologic descriptions and geophysical logs typically offer finer spatial resolution, and therefore more potential information about site-scale heterogeneity, than other site characterization-data; . . . . In this study, a technique for generating a heterogeneous, three-dimensional hydraulic conductivity field from sediment lithologic descriptions is presented. The approach involves creating a three-dimensional, finescale representation of mud (silt+clay) percentage using a "stratified' interpolation algorithm. Mud percentage is then translated into horizontal and vertical conductivity using direct correlations derived from measured data and inverse groundwater flow modeling. Last&, the fine-scale conductivity fields are averaged to create a 1 coarser grid for use in groundwater flow and transp or t modeling. The approach is demonstrated using a finite-element groundwater flow model of a Savannah River Site solid radioactive and hazardous waste burial ground. Hydrostratigraphic units in the area consist of fluvial, deltaic, and shallow marine sand, mud and calcareous sediment that exhibit abypt facies changes over short distances. For this application, the technique improves estimates of large-scale flow patterns and dispersive transport. The conductivity fields mimic actual lithologic data providing a more realistic picture of subsurface heterogeneity. Field-observed, preferential pathways for contaminant migration are replicated in the simulations without the need to artificially create zones of high conductivity. _INTRODUCTION Groundwater flow and especially contaminant transport can be affected by many physical, chemical and microbiological factors. Among these, large-scale (> 1 m) variability in hydraulic conductivity is the dominant influence on field-scale Characterizing hydrogeologic heterogeneity using lithologic data groundwater flow patterns and dispersive transport, with other factors introducing second-order effects (Brusseau, 1994). Incorporating hydraulic conductivity heterogeneity into flow and transport models is paramount to accurate simulations in most situations, particularly for contaminant migration. Sediment lithologic descriptions and down-hole geophysical logs typically offer finer spatial resolution, and therefore more potential information about heterogeneity, than most other site characterization data. However, these data do not provide direct information about hydrogeologic properties and can be viewed as "soft" hydrogeologic information. The emergence of powerful computing resources has enabled the routine use of spatiallydense, three-dimensional (3D) "soft" hydrogeologic data in 3D flow and transport models. Methods for supplementing "hard' hydrogeologic data (e.g., pump tests, slug tests, laboratory conductivity data), conventionally used to develop flow and transport models, with spatially-dense, 3D "soft" hydrogeologic data are the focus of this study. By utilizing widely available and 2 relatively inexpensive core data, these methods should reduce model uncertainty and the amount of field-scale hydrogeologic data required. Two key elements of these methods are 1) an algorithm for generating a complete 3D representation from discrete, spatiallyscattered field data, and 2) a method for translating "soft" data into hydraulic conductivity, and other necessary hydrogeologic parameters. Most articles in the literature focus on only one of these components. An exception is Lahm et al. (1995) who present a general procedure for estimating hydraulic conductivity from porosity (derived from geophysical logs) using a direct correlation, and generating a 3D conductivity field using a Bayesian updating/Cholesky decomposition stochastic approach. Below we review literghe falling in the two categories above, emphasizing studies most relevant to the geology of the US. Department of Energy's (DOE) Savannah River Site (SnS) (Fig. 1). Interpolation Algorithms Several methods have been proposed for generating a complete 3D representation Characterizing hydrogeologic heterogeneity using lithologic data from discrete, spatially-scattered field data including, deterministic interpolation, kriging, indicator geostatistics, turning bands methods, inverse and inverse conditional analysis, fractal-based approaches, and sedimentological models. All of these approaches have been used to successfully interpolate data. Each has inherent strengths and weaknesses, which should be considered when selecting an algorithm for a particular application. Deterministic methods such as inverse distance weighting and spline interpolation are conceptually simple, straightforward to implement, and generally yield good 3D representations of data. However, these methods typically ignore known or inferred information about spatial correlation, which can be used to improve the fit. Also, the uncertainty in interpolated values is difficult to estimate. For these reasons, investigators have pursued geostatistical techniques over the last two decades. Kriging (e.g., Journel and Huijbregts, 1978; Isaaks and Srivastava, 1989) originated from the mining industry as a geostatistical interpolation method that utilizes 3 information on the covariance of variables in space, in addition to data values. Kriging offers an improvement over simple interpolation schemes in many instances and has been widely employed over the last decade or so. The method also provides the uncertainty associated with each interpolated value, which may be important information depending on the application. Although an improved fit can be achieved, the choices of variogram and underlying error distribution are subjective and can lead to a wide range of possible solutions. Also, traditional kriging assumes that the covariance field is stationary or weakly stationary (only dependent upon distance between points), which may not be an appropriate assumption. In this case, additional effort is required to transform the original variable into a stationary one. Like other methods, kriging often perfo-poorly when there is a paucity of hard data. This is especially true for kriging methods that incorporate anisotropy (covariance dependent upon direction) to account for layering of sediments, for example. . Journel(l986) and Rubin and Journel (1991) attempted to overcome data Characterizing hydrogeologic heterogeneity using litholo& data deficiencies by utilizing "soft" information, and developed the binary indicator method with kriging to generate 3D fields. Phillips and Wilson (1989) and Brannan and Haselow (1993) also developed methods based on geometrical description of the geological media and statistical functions to generate fields. These turning bands methods provide additional refinement to the kriging approach, but can be labor intensive. Some natural processes are known to be fractal, and this observation has motivated recent investigations of fractal-based algorithms to generate 3D representations of geological systems. Hewett (1986) and Hewett and Behrens (1988) used fractal concepts to describe 3D geological systems with measurements of porosity from geophysical logs. Molz and Boman (1993) later applied the procedure to a hydraulic conductivity data base. Boman et al. (1995) concluded that while fractal-based methods are promising, their application requires further understanding and refinement. Some alternative methods generate 3D fields using knowledge of the physical processes of sediment deposition. For 4 example, Webb (1994) developed a simulator to mimic a braided stream depositional environment and generate 3D fields. These methods show promise, but will require signuficantly more research to encompass the complex and varied depositional processes that occur during the formation of a hydrogeological unit. Dagan (1985) proposed a method in which the 3D field is determined with the aid of an inverse solution of the groundwater flow equation. The method is termed a "conditional" or "unconditional" inverse simulation, depending upon whether measured data are honored in the interpolated field. Another inverse approach that has been proposed to generate 3D realizations is to solve the inverse problem for both the flow and transport equation to estimate the hydrogeologic properfies of an aquifer (e.g., Sun and Yeh; 1990; Wagner and Gorelick, 1987). Hyndman et al. (1994) proposed refining this method using seismic data. However, solving the flow and transport equations simultaneously is nontrivial and requires more development. Characterizing hydrogeologic heterogeneity using lithologic data Additional research is needed to evaluate and understand 3D interpolation algorithms in a practical sense, as Boman et al. (1995) have done in evaluating inverse distance weighting, kriging and fractal-based methods for generating a 3D conductivity field. Their evaluation indicates that a simple deterministic approach involving successive inverse distance weighting interpolation by horizontal layer is the preferred approach, based on its relative simplicity and the absence of improved tracer test modeling results using the other methods. Based on the conclusion of Boman et al. (1995), a relatively simple deterministic interpolation method was chosen for this study. Translation of "Sof t f Data into Hydrogeologic Parameters Methods employed in the petroleum industry for relating intrinsic permeability to borehole geophysical data are described in Lake and Carroll (1986) and Flint and Bryant (1993): Lahm et al. (1995) recently reviewed several efforts to relate permeability to borehole geophysical data including Timur (1968), Croft (1971), Keys and MacCary (1971), Kelly (1977), Biella et al. (1983), 5 Doveton (1986), Wendt et al. (1986), Urish (1981), and Jorgensen (1989). They also evaluated three methods of estimating conductivity primarily from porosity, as applied to the Milk River aquifer in Alberta, Canada. Lahm et al. (1995) found that a direct correlation between conductivity and porosity produced a lower residual between the measured and predicted datasets than the mefhods of Urish (1981) and Jorgensen (1989). Many predictive methods for estimating hydraulic conductivity from grain-size data have been proposed over the last several decades Freeze and Cherry, 1979; Domenico and Schwartz, 1990). These methods are most reliable when applied to aquifer zones consisting of clean sand and gravel. Riha (1993) evaluated ten predictive equations as applied to Coastal Plain sediment bjeneath the SRS and found that none performed adequately. The methods evaluated were those of Hazen (1892), Kozeny-Carmen (in Carmen, 1937), Krumbein and Monk (1943), Harleman et al. (1963), Masch and Denny (1966), Beard and Weyl(1973), Puckett et al. (1985), Ahuja et al. (1989), Franzmeier (1991) 8 / Characterizing hydrogeologic heterogeneity using lithologic data and Jabro (1992). Riha (1993) developed a predictive equation from a small SRS data set, but that correlation has not been tested beyond the data used to develop the correlation. Kegley et al. (1994) performed over 4000 minipermeameter tests of Tertiary-aged Coastal Plain sediments at the SRS and correlated the geometric mean of measured permeability for each stratigraphic unit within each well to mud fraction estimated by binocular microscope examination. The correlation coefficient (r) of the expression is 0.92. Present Study For this study, borehole geophysical data are combined with information from microscopic core examinations to yield qualitative grain-size distribution information in the form of estimated percent gravel, sand, and mud. A general approach is presented for generating a heterogeneous, 3D hydraulic conductivity field from estimated mud fraction (silt + clay sized material). The mud fraction data are first interpolated onto a 3D Cartesian grid using the minimum-tension spline algorithm 6 implemented in Earthvision@ version 2.0 (Dynamics Graphics, Inc., Alameda, CA) with minimal “vertical influence” selected. This interpolation approach is analogous to the successive inverse distance weighting interpolation by horizontal layer scheme recommended by Boman et al. (1995), and is easily applied using Earthvision@. The mud fraction grid is translated into horizontal and vertical conductivity fields using direct . correlations of conductivity to mud fraction. The correlations are developed initially from laboratory data and adjusted through inverse groundwater flow modeling. Selection of a direct correlation approach is motivated and supported by Kegley et al. (1994) and Lahm et al. (1995). The fine-scale conductivity fields are then transferred to a coarser grid for use in groundwater flow modeling through arithmetic and harmonic Weraging. This method is demonstrated in a finiteelement groundwater flow model for a radioactive and hazardous waste burial ground at the Savannah River Site. SITE DESCRIPTION AND HYDROGEOLOGY The study area is located in the Upper Atlantic Coastal Plain physiographic Characterizing hydrogeologic heterogeneity using lithologic data province in southwestern South Carolina at the U. S. Department of Energy Savannah River Site (SRS) (Fig. 1). The SRS occupies approximately 800 km2 and was set aside in 1950 as a controlled area for the production of nuclear material for national defense, and specialized nuclear materials for medical and industrial purposes. The study area is located near the center of the SRS within the 40 km2 General Separations Area (GSA) (Fig. 1). This area is undergoing intense hydrogeological characterization for RCRA and CERCLA investigations. The focus of this study is the original solid radioactive and hazardous waste burial ground in the GSA. The burial ground received solid radioactive and hazardous waste from approximately 1952 to 1972. Primary groundwater contaminants include tritium, trichloroethylene (TCE) and tetrachloroethylene (PCE). Groundwater plumes emanating from the burial ground are migrating to the south toward Fourmile Branch (Fig. 1). The SRS is underlain by a southeastdipping wedge of unconsolidated and consolidated sediment consisting of sand, mud, limestone, and gravel, which range in age from early Late Cretaceous (Cenomanian) to Holocene (Fallaw and Price, 1995). This study involves the Tertiary-aged sediment, principally the Eocene to Miocene sequence (Fig. 2). The Gordon aquifer, the Gordon confining unit and the Upper Three Runs aquifer are the hydrostratigraphic units of interest (Fig. 2). The Gordon aquifer consists of unconsolidated sand and clayey sand of the Congaree Formation and, where present, the sandy parts of the underlying Williamsburg Formation. The sand is yellowish to grayish orange in color and is subrounded to wellrounded, moderately to poorly sorted, and medium to coarse-grained. The Gordon confining unit separates the Gordon and Upper Three Runs aquifers and is lithostratigraphically equivalent to-the Warley Hill Formation. This unit is composed of fine-grained, silty and clayey sand, sandy clay and clay. The clay is stiff to hard and sometimes fissile. Zones of silicacemented sand and clay, and calcareous sediment are common. The Upper Three Runs aquifer includes